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OPTIMAL EXPERIMENT DESIGN IN A FILTERING CONTEXT 
WITH APPLICATION TO SAMPLED NETWORK DATA 

By Harsh Singhal and George Michailidis 1 

University of Michigan, Ann Arbor 

We examine the problem of optimal design in the context of fil- 
tering multiple random walks. Specifically, we define the steady state 
E-optimal design criterion and show that the underlying optimiza- 
tion problem leads to a second order cone program. The developed 
methodology is applied to tracking network flow volumes using sam- 
pled data, where the design variable corresponds to controlling the 
sampling rate. The optimal design is numerically compared to a my- 
opic and a naive strategy. Finally, we relate our work to the general 
problem of steady state optimal design for state space models. 

1. Introduction. Consider a wide area computer network such as the 
one depicted in Figure 1. A flow is defined as all traffic with common ori- 
gin and destination nodes. Monitoring flow volumes plays an important role 
in network management tasks, such as capacity planning by tracking de- 
mands and forecasting traffic, identifying failures together with their causes 
and impact, detecting malicious activity and configuring routing protocols 
[Barford et al. (2002), Soule et al. (2005)]. These flow volumes have been 
observed to exhibit complicated structure, as seen in Figure 2. For example, 
the highly aggregated flows usually have diurnal patterns [Figure 2(a)], while 
lighter flows can be extremely noisy [Figure 2(b)]. Network traffic is carried 
on packets that can be observed (and sampled) at router interfaces, hence- 
forth called observation points. However, during the measurement process, 
sampling is employed due to high flow volumes and resource constraints at 
routers. 

It is increasingly common for such measurement infrastructure to be de- 
ployed in computer networks [Duffield (2004)]. Each packet from the ag- 
gregate flow at an observation point is sampled independently with a cer- 
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Fig. 1. Geant Network: (a) geographic view (www.geant.net) and (b) the corresponding 
logical topology. 



tain probability (sampling rate) [Duffield, Lund and Thorup (2004)]. Typi- 
cal sampling rates range between 0.001-0.01. For every packet sampled, its 
header information is recorded which allows one to reconstruct objects of 
interest, such as volumes of flows with a particular source and destination 
traversing the network. An important issue is how to select (design) the 
sampling rates across the network subject to resource constraints, in order 
to collect the maximum amount of information on the underlying source- 
destination flows. Obviously low sampling rates result in large sampling 
noise. One way of achieving lower estimation error with the same sampling 
rate is through filtering; that is, combining the present measurement with 
past measurements to track the time-series of flow volumes. In designing a 




Fig. 2. Flow volumes: (a) all flows and (b) one of the lighter flows. 
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sampling scheme for this situation one needs to take into account measure- 
ment noise and process noise (innovation noise). 

While modeling the dynamics of flow volumes is a challenging task in 
itself [Park and Willinger (2000)], we use a simple random walk model for 
this purpose. This is a robust enough model to be useful in a large range 
of applications and leads to scalable filters. We consider the problem of 
minimizing the (running) estimation error through optimal design of mea- 
surement scheme in the filtering context. In this paper we take an optimal 
design of experiment approach to the above problem and demonstrate its 
application to computer network monitoring using sampled data. 

The related research on optimal design has focused on one of the following 
scenarios. There is a large body of work on optimal input design for dynam- 
ical systems [Goodwin (1977), Titterington (1980)]. There the focus is on 
parameter estimation (system identification) rather than filtering, as in this 
paper. Another related area is sequential design for nonlinear systems [Gau- 
tier and Pronzato (1998), Ford, Titterington and Kitsos (1989)], where the 
optimal design depends on values of unknown parameters. While there are 
some commonalities, the design problem in a filtering context is unique in 
that the design at any time affects not just the current estimation error but 
also future ones. The problem of optimal sensor placement in control sys- 
tem literature looks at an equivalent problem [Arbel (1982), Chmielewski, 
Palmer and Manousiouthakis (2002)]. However, the formulation is not in 
terms of information matrices and the special case of random walks has not 
been analyzed to our knowledge. More details are provided in Section 4. 

The remainder of the paper is organized as follows: in Section 2 we formu- 
late and investigate the idealized problem of optimal design in the context 
of filtering for multiple random walks. In Section 3 we study its application 
to tracking flow volumes using sampled data. We end with discussion of a 
possible generalization and some comments in Section 4. 

2. Optimal design for multiple random walks. Let us first briefly review 
the concept of E-optimality from classical design-of-experiment literature 
for a simple setting. Assume we have independent observations 



for i = 1,2, . . . ,n r . The natural estimate for X{ is = for all i. It is 
standard to assume that the inverse variance of observation noise is roughly 
proportional to design variables. The inverse variance, m^, can be thought 
of as the information collected on parameter X{. Specifically, we assume that 
the relation between an n r x 1 information vector m and an n Q x 1 vector 
of design variables £ is 





(2.2) 



m = J£. 
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For example, suppose there is a library of measurements z±, . . . ,z no , each 
of which is independently distributed as Zj ~ N(x^,afl), where is a 
subset of elements of x. Let £j be equal to (or proportional to) the number 
of independent measurement of type i (replications of Zi) collected during 
the experiment. Then, the weighted least squares estimate y of x can be 
shown to have distribution given by (2.1) and (2.2). The matrix J depends 
on the the membership of subsets xui and variances of (assumed known), 
for £ = 1, ... , n Q . 

We assume that the design variables are constrained to be positive and, 
in addition, satisfy n v linear inequality constraints. These can be written 
as R^ < b, where R is an n v x n Q matrix and b is n v x 1 vector. We think 
of this type of constraint as a budgetary one, that specifies upper limits on 
weighted sums of the design variables. 

Now the E-optimal design problem is given by 

are max min rrij . 

The objective function, minimi, is the minimum information over all 
flows. Note that this corresponds to minimizing the maximum mean squared 
error (MSE) since 1/rrii is the MSE in the estimate of Xj. Using maximum 
MSE as the objective function corresponds to aiming for the best possible 
worst case performance. 

As an example consider the situation where m\ = 40£i + 10^2 and m-2 = 
10£i +40^2- Further assume the constraint 
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Fig. 3. Contours of the objective function for E-optimal design. 
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Figure 3 shows the contours of the objective function, that is, minjTOj. 
The region below the thick line is the constraint space. As usual, the optimal 
solution corresponds to the point where the contour of the objective function 
is "tangent" to the boundary of the constraint space. It is clear that the 
optimal design would be £1 = £2 = 0.50, which is also reasonable from the 
symmetry of the setup. 

We now extend the above criteria to the optimal design for random walks 
in a steady state. Consider a collection of independent random walks 

Xi(t) =Xi(t-l)+£i(t), 

for i = 1, . . . ,n r and t = 1, 2, The term £i(t) is called the innovation noise 

and we assume that Var(£j(i)) = erf, which is referred to as the innovation 
variance. Further, suppose we have noisy observations 

Ui(t) =Xi(t) +rji(t). 

Let Vax(rji(t)) = 1/mj. As before, we assume the relation between observed 
information and design variables to be m = J£, with n r x n Q matrix J as- 
sumed known. 

The estimates of interest in this case are the ones obtained through fil- 
tering 

x i (t)=E[x i (t)\y i (t),y i (t-l),...]. 

Let Si(t) = Vax(xi(t)\yi(t),yi(t - 1),...). Further, let fhi = limj^oo l/s;(t) 
when it exists. We will refer to this as the steady state information. When 
the innovation and measurement noise, £i(t) and r]i(t) respectively, are Gaus- 
sian, the optimal filter corresponds to a Kalman filter and in this case the 
steady state always exists [Harvey (1990)]. For the remainder of the paper 
we will assume that £i{t) and rji(t) are independent mean Gaussian random 
variables. If Si(t\t — 1) = Var(xi(t)\yi(t — l),yi(t — 2), . . .), then the Kalman 
filter update equations give us 

(2.3) Si (t\t - 1) = Si (t - I) + af 
and 

(2.4) Siit)- 1 =s i (t\t-l)~ 1 +mi 
(2-5) =( u \w. +4\ + m i- 



Si(t-l) 

Note that, given Sj(0), of and rrii, one can calculate Si(t) at any time t by 
iterating the above equations. Further, the choice of m; impacts Si(t) not 
just for a specific t but for all t. Thus, 

m = ( ) + rrn 

rrn 
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Fig. 4. Contours of the objective function for steady state E-optimal design. 



or 

2-2 2 n 

a i m i — a i miUii — rrii = 0. 

Hence, 



We define the steady state E-optimal design problem as 

argmaxmin rrij. 

R£<b i 

As an example consider the same setting as above, with mi = 40£i + 10^2 
and m,2 = 10^i + 40^2- Further, let the innovation noise be characterized by 
o"i = 0.1 and 02 = 0.2. As before, we assume the design constraint 

6 + 6<i. 

Notice in Figure 4 that even though there is symmetry in the measured 
information, the first random walk is smoother than the second one and, 
hence, less measurement resources need to be allocated to it. 



2.1. Optimization for the steady state E-optimal design. We establish 
next the main technical result of the paper, that the steady state E-optimal 
design problem is a second order cone program. First, we introduce a new 
variable 9 as the lower bound for the steady state information over all flows. 
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To solve the steady state E-optimal design problem, we have to maximize 9 
subject to 



( 2 -6) Hr-j > 



m i a i + \/ m i °f + 4m jcri 

~R 

for i = 1, . . . , n r and 

R£<b. 

Equation (2.6) can be equivalently written as 

(2.7) e 2 <m i (e + X 1 

which is a hyperbolic constraint [Lobo et al. (1998)]. Thus, this problem can 
be cast as a second order cone program (see the Appendix for a review of 
second order cone programs and the representation of the above optimization 
in canonical form). Such optimization programs can be solved efficiently 
through interior point methods [Boyd and Vandenberghe (2004)], software 
implementations of which are commonly available [Benson and Ye (2008), 
Grant and Boyd (2009)]. 

2.2. Myopic approach. In the following, we present a greedy alterna- 
tive to the steady state optimal design. As before, assume yi(t) = Xi{t) + 
r?j(i). Further, we assume that Var(r/j(t)) = l/mj(t); that is, we allow for 
time varying design variables £(t) with m(t) = J£(t). As before, Sj(t) = 
Var(xi(t)\yi(t) , yi(t — 1), . . .). Define the information at time t to be given by 
rhi(t) = l/si(t). Note that fhi(t) is a function of — 1), 

The myopic E-optimal design at time t is defined as 

are max minm,(t). 

mt)<b i 

Note that the objective function only involves rhi(t), that is, the infor- 
mation at time t. However, the choice of impacts not just rhi(t) but 
also fhi(t + 1), rhi(t + 2), . . . , due to the iterative nature of Kalman filtering. 
Since it ignores this "long term impact," we refer to this scheme as myopic. 
Equation (2.5) implies that 

mi(t) = Si (t\t - I)" 1 + J£(t). 

As before, a new variable 9 can be introduced to lower bound rhi(t) which 
gives a new set of constraints 

Si (t\t-iy l + Jt(t)>o, 

in addition to the original constraint R£(t) < b. Now the objective is to 
maximize 9 with the optimization variables being 9 and £(t). Since both 
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Fig. 5. Schematic representation of an observation, observation points and flows. 

the objective function and the constraints are linear in and 8, the above 
optimization is a linear program. Not surprisingly, the myopic optimal design 
is a much easier problem than the steady-state optimal design even in more 
general settings as noted in Section 4. Note that since the sampling rates are 
allowed to vary with time, the myopic optimal design may have an objective 
function larger than the steady state optimal case. However, as the objective 
of optimization is to maximize present information with no regard to impact 
on future information, such a scheme can not be guaranteed to perform well 
in the long run. 

3. Application to tracking flow volumes. The ideas developed above can 
be used for designing the sampling rate in a computer network for tracking 
flow volumes. As mentioned in the introduction, we will use the random 
walk model for flow volumes due to its simplicity and robustness. 

Suppose there are n r origin-destination flows in a network. Let Xi(t) be the 
volume of the ith flow in the t-the time interval, for % = 1, . . . , n r . These flow 
volumes are tracked using sampled data which are noisy. Recall that flows 
are sampled at router interfaces, which we refer to as observation points. In 
the past, a systematic sampling scheme was the dominant technology, but 
truly random sampling technologies have recently become available and are 
commonly deployed [Driffield, Lund and Thorup (2004)]. All flows travers- 
ing an observation point (router interface) experience the same sampling 
rate. Each incoming edge at a node in Figure 1(b) is an interface of the 
corresponding router. Each router typically has multiple interfaces and each 
flow may traverse multiple observation points due to multi-hop paths and 
multi-path routing. 

Suppose there are n a observation points on the network where sampled 
data on flows can be collected. Further, assume that sampling rates of 
£ = (£i, . . . ,£n )' are used at observation points 1, . . . ,n a , respectively. Any 
given observation point k € {1, . . . ,n } generates estimates for elements 
of x(t), that is, the number of flows that go through that node. Thus, a 
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total of n g = Ylk=i 9 k measurements are available in each time interval t, 
say, Z\(t), . . . ,z ng (t), which need to be optimally combined to get the re- 
quired estimates. Consider the router in Figure 5. Assume that k(i) is the 
observation point at which the ith measurement is collected, i = 1, . . . ,n g , 
and the corresponding flow. Thus, k(-) : {1, . . . ,n g } — > {1, . . . ,n } and 
/(•) : {1, . . . , n g } — > {1, . . . , n r }. Further, let 

E[zi{t)\x l{i) {t)} = x l(j) {t) 

and for the moment assume 

(3-1) Cov(zi(t)\x l{i) (t)) = Hi(i)/£ k (i), 

where /i, = E[xi(t)]. The exact sampling mechanism and approximation in- 
volved in the above relation are described in Section 3.2. Thus, in vector 
notation we get 

(3.2) E[z(t)\x(t)]= Lx(t), 

where L is a n g x n r matrix with Lij = 1 only if = j (i.e., ith measurement 
corresponds to j flow) and otherwise and 

(3.3) Cov(z(t)\x(t)) = D, 

where D is a n g x n g diagonal matrix, with [D]u = fJ-i(i) / S,k(i) ■ Using (3.1), 
the inverse of D is given by D^ 1 = X^fc^fc^fc) where ^fc, k = 1, . . . ,n , are 
rig x rig diagonal matrices with their zth element given by 

Let y(t) be the general least squares estimate of x(t), under equation (3.2) 
and (3.3). Thus, 

(3.5) Cov{y(t)\x(t)) = (L'D^L)- 1 

(3.6) = (^(L'* fe L)£ fc 

^ k 

From the definition of L, it follows that the jth elements of any two 
columns of L cannot be nonzero simultaneously. Thus, the matrix in (3.5) 
is diagonal. Further, 

DiagtCov^t)!^)))- 1 =?n= J£, 

where 

[J] ik = L'^ k L. ti . 
We will refer to the above as the linear model. 
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Sampling is employed in network flow measurements because measure- 
ment resources like CPU time and available storage are limited. Typically, 
all observation points (router interfaces) belonging to a particular router 
share these resources. We assume that the sampling rates are constrained 
to lie in a convex polygon R£f < b. This includes the case where the sum of 
sampling rates on the interfaces of a router is bounded above by the budget 
for that router. We will focus on this constraint for the rest of the paper. In 
this case, the constraints are given as one linear inequality for each router. 

For the available data, we set up the performance evaluation as follows. 
We use the Geant network topology, which has n v = 23 nodes (routers) and 
37 x 2 bidirectional edges. The available data [Uhlig et al. (2006)] correspond 
to flow volumes over time. Each time interval is equal to 15 minutes. The 
original data set spans 4 months, but we focus on the first 200 time intervals 
to avoid severe non-stationarities inherent in an evolving network. Further, 
we focus on the top 25% of measured flows by volumes since one is typi- 
cally interested in tracking heavy flows. This corresponds to n r = 76 flows. 
We assume that sampled data can be collected at each incoming edge of a 
router and, thus, we have n Q = 37 x 2 observation points. We assume that 
these flows are routed through minimum distance paths, which is a common 
routing mechanism in wide area networks [Peterson and Davie (2003)]. This 
leads to n g = 163 and the routing information gives us the mapping l(-) and 
hence the matrix L. Matrix L is 163 x 76, matrices D, ^f\, . . . , ^74 are all 
163 x 163. We assume that the sum of sampling rates on all interfaces of a 
router is bounded above by 0.01, that is, bi = 0.01 and Rij = 1 if observa- 
tion point j is an interface of router i and otherwise. Thus, matrix R is 
23 x 74. Finally, we estimate the af and \%i parameters associated with the 
flow volume processes, and assume they are available for filtering purposes 
and measurement design. As we have argued, both the steady state optimal 
and myopic design problems are standard optimization programs and once 
they are written as such, any standard optimization package [Benson and 
Ye (2008), Grant and Boyd (2009)] can be used to solve them numerically. 

For the purpose of comparison, we also define a naive sampling scheme as 
follows. For any given router, an equal sampling rate is allocated to every 
interface that carries any of the 76 flows of interest. This allocation is done so 
as to make the corresponding budget constraint tight. For example, suppose 
the ith router has 5 interfaces, but only 4 of them are traversed by one 
of the 76 flows of interest. In this case, each of the latter 4 interfaces will 
be allocated a sampling rate of 6^/4, while the remaining interface will be 
allocated a sampling rate of 0. 

3.1. Performance of various sampling schemes for the linear model. Fig- 
ure 6(a) shows the value of the maximum MSE as a function of time. Note 
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that as information accumulates over time, we obtain an improvement in per- 
formance under all three sampling mechanisms, myopic, naive and steady 
state optimal. Here performance is measured as the maximum of Si(t) over 
all flows, calculated using equations (2.3) and (2.5). Surprisingly, both the 
myopic and steady state optimal sampling mechanisms perform equally well 
in the steady state and achieve a 42% improvement over the naive sampling 
in the steady state. Figure 6(b) shows that the myopic optimal sampling 
rates at all observation points reach a steady state. Figure 7 shows the value 
of steady state sampling rates at various router interfaces in the network 
topology. Even though the myopic scheme has the flexibility of time varying 
sampling rates, if the sampling rates do reach a steady state, its performance 
can clearly be no better than the steady-state optimal scheme. However, as 
Figure 6 shows, in this case, the additional flexibility permits the myopic 
scheme to reach steady state performance faster than the steady state opti- 
mal one. 



3.2. Departures from the linear model: performance with geant data. A 
more detailed model for flow volumes and sampled measurements would 
have to include significant departures from the linear model assumed above. 
First, the true flow volumes clearly have more structure than independent 
random walks, as seen in Figure 2. In applying the above ideas to the Geant 
data, we will investigate their robustness to the independent random walk 
assumption. 

A more serious departure is the following. Suppose that a flow with volume 
X in a certain time interval is sampled at a rate £. If the number of sampled 
packets is N, then the usual (approximate maximum likelihood) estimate of 
flow volume is Z = N/£. The variance of measurement noise can be shown 
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Fig. 6. Performance of various sampling schemes [panel (a)/ and sampling rates at var- 
ious interfaces under a myopic scheme [panel (b) ]. 



12 



H. SINGHAL AND G. MICHAILIDIS 



0.01 -, 



0.008 - 



0.006 - 




Fig. 7. Spatial view of steady state optimal sampling rates. 



to be Var(Z|A) ~ X/£ [Duffield, Lund and Thorup (2002)]. Thus, m in (3.1) 
is actually equal to the unknown Xi(t). 

The observation above implies that in applying the presented techniques 
to sampled network data, one would have to rely on an approximate model 
for measurements Zi (t) . We will follow an approach similar to batch sequen- 
tial design [Gautier and Pronzato (1998)]. Assume that the sampling rates 
are to be held constant for a batch of contiguous time intervals. At the be- 
ginning of each batch, we use the most recent estimate Xi(t — 1) in place 
of fa in (3.4) for sampling rate design. For filtering purposes, we employ a 
Kalman filter with Xi(t — 1) in place of jjn in equation (3.4) at each time t. 
We replace the budget constraint inequalities R£ < b with the corresponding 
equalities R£ = b to force full utilization of available resources. For routers 
that are traversed by at least one of the 76 flows of interest, we introduce ad- 
ditional equality constraints as follows. Design variable for an interface k 
not traversed by one of the 76 flows of interest is constrained to be identically 
0. Figure 8 shows the performance of different sampling schemes averaged 
over 200 realizations of sampled data. The sampled data emulate the exact 
sampling mechanism described above (with respective sampling rates) with 
the Geant data treated as the underlying (unobserved) flow volumes. 

Sampling rates were adjusted only at the beginning of a 40 time period 
block and were held constant over each block. In the first block, the sam- 
pling rates were forced to be the same as the naive scheme irrespective of the 
sampling mechanism under study. Notice that for low values of the objec- 
tive function (maximum mean squared error) the myopic and steady state 
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Fig. 8. Performance of various sampling schemes using batch sequential design with flow 
volumes from the Geant data. 

allocations perform better than the naive allocation. On the other hand, 
when the maximum mean squared error spikes, the naive allocation per- 
forms better, indicating robustness to model departures. The median (over 
time periods 41 to 200) of maximum MSE for myopic, naive and steady-state 
optimal sampling is 5.46 x 10 9 , 7.49 x 10 9 and 6.14 x 10 9 , respectively. Thus, 
the myopic scheme performs better than the steady state optimal scheme, 
which in turn performs better than the naive scheme. 

Finally, we look at the performance of the myopic allocation when the 
above scheme is employed with a block size of just one time interval; that is, 
sampling rates were adjusted at the beginning of each time period using the 
myopic scheme. The results are displayed in Figure 9. As before, the current 
estimate of flow volumes is used in place of /Xj in equation (3.4) for both 
filtering and myopic sampling scheme design. The myopic sampling scheme 
can be seen to perform better than the naive version in most time periods. 
The median (over time periods 1 to 200) of maximum MSE is 4.50 x 10 9 
and 7.44 x 10 9 for myopic and naive sampling, respectively. 

4. Discussion and future work. The specification of the steady state op- 
timal design problem can be easily generalized to linear dynamical systems. 
Such systems are described by a pair of equations [Harvey (1990)]. The state 
transition equation can be written as 

x(t) = Cx(t-l)+w(t), 

where Cov (w(t)) = W. The observation equation can be written as 

y(t)=Lx(t) + e(t). 
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Fig. 9. Performance of fully time varying myopic and naive sampling mechanism with 
flow volumes from Geant Data. 

Assume Cov(e(t)) = , where is a linear function and £ is the 

value of design variables. 

For the above dynamical system a Kalman filter can be used to iteratively 
compute E[x(t)\y(t),y(t — 1), . . .]. Let the steady state estimation error co- 
variance be £ = M~ l (assuming the system is observable [Harvey (1990)]). 
Then, M satisfies the Algebraic Riccati equation: 

(4.1) M={CM~ 1 C' + W)~ l +L'V(£)L. 

Such equations have no analytic solution in general. 

The steady state optimal design problem can now be defined as 

argmax f(M), 
R£<b 

where /(■) is an appropriate scalarization of the information matrix [Fedorov 
and Hackl (1997)]. An interesting open problem is to solve the above opti- 
mization efficiently in the absence of an analytic solution to (4.1). The sensor 
placement problem in control system literature [Arbel (1982)] is equivalent, 
though not identical. The Newton- type algorithm proposed in [Arbel (1982)] 
for this problem requires the solution of the Algebraic Riccati equation at 
each iteration of the algorithm. It would be desirable to develop more effi- 
cient algorithms. 

In summary, we have shown that steady state E-optimal design for ran- 
dom walks is a second order cone program. We have illustrated numerically 
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that the performance of the Kalman filter can be significantly improved by 
incorporating an optimal experimental design. The linear state space model 
is of general interest and one would like to investigate the steady state opti- 
mal design problem described above. Finally, from a practical point of view, 
it would be useful to extend these ideas to nonlinear filtering. 

APPENDIX: OPTIMIZATION REVIEW 

In this section we summarize the concepts of second order cone programs 
and hyperbolic constraints from Lobo et al. (1998). We also present the 
steady state optimal design problem in the canonical form. 

A second order cone program is defined as 



minimize fx 

subject to \\PiX + qA\ < r[x + Sj 



1 N. 



Here, x G W 1 is the optimization variable, and the problem parameters are 
/ G M n , Pi 6 M n ' xn , qi G M ni , rj G R n and Si G R. The norm in the constraints 
is the standard Euclidean norm. A second order cone program is a standard 
convex program and algorithms to numerically solve it are well studied and 
implemented in computational software. 
A constraint of the form 

w < xy, x > 0, y > 

is called hyperbolic. The above can be shown to be equivalent to 

2w 
x-y 



<x + y. 



Using the above representations, we can write the steady state optimal de- 
sign problem as a canonical second order cone program as follows. Equation 
(2.7) can be equivalently written as 



m; 



1 



07 



< nii + 9 + 



1 



07 



Thus, the steady state optimal problem is a second order cone program 
with N = n r + n v , a/ = (0,&,...,£nj, /' = (-1, 0, . . . , 0), 



1 





Ji,l 



Ji. 



l,...,n r and Pi = 0, for i = n r + 1, 



,n r + n v , q\ 



(0,-1/ff?), for 



for i 

i = l,...,n r and % = 0, for i = n T + 1, . . .,n r + n v , r[ = (1, J it \,.. . , Jj, no ), for 
i = 1, . ..,n r and r- = (0, -Ri- Urj i, . . ■ ,-Ri-n r ,n ), for i = n r + 1, . . .,n r + n v , 

h 



Si = l/<7~, for i = 1, . . . , n r and Sj 



., for i = n r + 1, . . . ,tv + n„. 
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